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ABSTRACT 


Results are obtained for ?-D hygroscopic diffusion into composite 
laminates with free edges through the use of a second order finite 
difference scheme. A computer program was developed for the transient 
nonlinear coupled hygrothermal diffusion into finite width laminates 
with varying surface conditions along two edges. The formulation 
permits the diffusion coefficient to be a function of temperature, 
moisture concentration, and fiber orientation. The moisture distribu- 
tions thus obtained are necessary for analysis of the moisture induced 
interlaminar stresses in the vicinity of free edges. It is shown that 
large diffusion times tend to eliminate the significance of stacking 
sequence, but gradients within and between layers are significant for 
all diffusion times. 
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INTRODUCTION 


Problems associated with combined thermal and mositure effects in 
polymeric matrix composite materials are currently receiving consider- 
able attention on the part of the materials research community. The 
experimental works of McKague et al [1] and Browning et al [2,3] have 
shown that there is a degrading effect on the shear, compression and 
transverse stress properties of resinous composites due to moisture 
absorption. These degrading effects are primarily due to changes in 
properties of the resinous matrix and the fiber-matrix interface. 

McKague et al also presented rather extensive experimental results 
showing that the rate of moisture absorption is temperature dependent 
with the rate of absorption increasing with increasing temperature. It 
was also shown in reference [3] that large moisture concentration 
gradients may exist near the surfaces of composite laminates. These 
gradients may be significant because of the possibility of large inter- 
laminar stresses near the free edge of composite laminates. 

Analytical treatments of moisture diffusion in composites have been 
presented by several authors. Shen and Spring [4] considered one- 
dimensional absorption and desorption of homogeneous materials and 
composites, and Whitney [5] considered three-dimensional moisture 
diffusion in laminated compo.sites. The emphasis in both of the above 
papers was on the total percent moisture weight gain as a function of 
time. Whitney also presented moisture profiles for two-dimensional 
diffusion into a laminate. These profiles were generated assumming the 
diffusion coefficient was an averaging of the laminates constituents. 
Since edge effects were included through a correction factor and were 
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not treated explicitly in either of these two later papers, the signi- 
ficance of laminate stacking sequence wasn't readily assessable. 

This paper presents numerical results for two-dimensional moisture 
diffussion into a symmetric finite width composite laminate. The 
results were obtained from a recently developed finite difference 
computer capability, denoted HYDIP, and the associated graphics program 
HY0IP6. This program was developed for the two-dimensional solution of 
the transient nonlinear coupled hygrothermal diffusion into a symmetric 
finite width composite laminate with varying surface conditions along 
two edges. The moisture distributions thus obtained are necessary for 
analysis of the interlaminar stresses which may be present in a com- 
posite laminate with free edges since these interlaminar stresses often 
initiate failure. 

2. PROBLEM FORMULATION 

The governing 3-D field equation for hygroscopic or thermal dif- 
fusion into a body is given by 

3 /u 3R\ .3 tu 3R\ . 3 (u 3R\ _ 3R /-i\ 

3x ' X 3x 3y l *y 3y ; 32 3Z j " 3t []) 

where K , K and K represent the diffusion coefficients in the x, y and 
x y z 

z directions, respectively, and R is the temperature or moisture con- 
centration at a point in the body. Restricting the problem to that of a 
finite width symmetric laminate wih no variation in the x direction, 
figure (1), only a quarter laminate need be analyzed. The governing 
field equation becomes: 

3 / w 3R \ , 3 / 1/ 3R x 3R 

3y ' y ay' 3z ' z az' at 


( 2 ) 
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with symmetry boundary conditions 


8R 

3Z 


= o- 22. 
2=0 8 > 


y- o 


= o 


(3) 


and variable surface conditions 


R = R(z) b ; R ■ R(Y) z= 


H‘ 


(4) 


More specifically for the case of hygroscopic diffusion, the dif- 
fusion coefficient 0 may be a function of temperature T(y,z), moisture 
concentration M(y,z) and lamina fiber orientation e, 

D = D(T(y,z) , M(y,z) , fl). 

The governing field equations, (2) - (4), can now be written as 


^ M + D ^ + 

ay ay y 3y 2 


3H n 3 2 M _ 9M 

3Z dZ U Z „,2 " 3t 

oZ 


(5) 


with symmetry boundary conditions 


3M 

9Z 


= o* M 

z-0 ^ 


= 0 


y=o 


end variable surface conditions 


( 6 ) 


M = M(z) ; M - M(y) Z=H - 


(7) 


3. SOLUTION SCHEME 

The solution of the diffusion equation (5) was performed by a 
second order accurate variable step size finite difference scheme [6]. 
First and second derivative central difference derivatives are employed 
in the interior region of the laminate while only first derivative 
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forward difference derivatives are necessary along the symmetry boundary 
conditions. A complete presentation of the basic difference equations 
used and the equivalent difference equation of equations (5) and (6) are 
presented in Appendix A. 

Applying equations (All) - (A13) at nodes throughout the laminate 
yields a set of linear simultaneous equations of the form 

* {P t> < 7 > 

where [C t ] is the coefficient matrix at time t, are the unknown 

moisture concentration at time t+1 and {P^} is the right hand side 
vector which is a function of the moisture concentration at time t. It 
should be realized that at t * 0 {P^.q) would be a function of the 
initial moisture distribution. Equation (7) is then successively solved 
over N time increments. Transient problems can similarly be solved 
by imposing transient surface conditions in an incremental fashion. 

An iterative simultaneous equation solution scheme [6], though more 
expensive to use for a single analysis, was chosen over an elimination 
or decomposition ilgorithm. This scheme was selected due to the in- 
cremental nature of the diffusion problem and because iterative schemes 
become cost efficient when successive solutions change only slightly 
from the previous results. 

A Gauss-Siadel simultaneous equation solver was modified to ef- 
ficiently solve equation (7) over N time intervals. This was primiarily 
achieved by eliminating the storage of zero coefficients and the multi- 
plication by zeros. The required storage locations and multiplications 
of the modified Gauss-Siedel routine then become linear functions of 
the number of unknowns. A more detailed discussion of this solution 
scheme is presented in Appendix B. 



4. CASE STUDIES 


As mentioned in the introduction, large gradients in mositure 
concentration may exist near the free surfaces of a laminate. In this 
study it will be assumed that these gradients are the result of the 
laminate surface being subjected to an environment different from the 
one associated with the moisture distribution initially existing in the 
laminate. This condition could similarily be induced by thermal gradients 
in the composite as well as a variety of combi nations of thermal and 
mositure distributions. 

The initial studies investigated the effect of stacking sequence on 
the moisture distribution within the boundary layer region, approxi- 
mately two lamina thicknesses along the laminates free edge. Material 
diffusion coefficients, [2], used in this study are shown in figure (2) 
while the laminates investigated are: [90 2 ,0 2 3 s , [90,0 3 ] s , [0,90 2 ,0] $ , 

[90,0 2 ,90] s , [0 2 ,90 2 ] s , [0,±45,90] s , [90,±45,0] s> 

Upon examination of the laminates studied it becomes obvious that 
primarily 0 N /90 N combi nation laminates were considered. These laminates 
were chosen because they produce the most dramatic results due to the 
relative magnitude of the diffusion coefficients. Though the studies 
were very selective, the program HYDIP will handle any number of layers 
each at differenv fiber orientation and each having different material 
properties. 

The basic model was a four layer symmetric laminate with 510 nodes 
on the model. Each layer, 0.005 in. thick and 1.0 in. wide, contains 
four nodes through the thickness with an additional node on the upper 
surface and 30 linearly spaced nodes across the width. Figure (1). The 
laminate was assumed to be initially dry and at a temperature of 440°K. 
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Tills temperature was chosen because the diffusion coefficient was much 
larger than that at room temperature as seen In Figure (2). An applied 
surface condition corresponding to 100 percent relative humidity was 
uniformly distributed along both free surfaces. 

The plotted results In figures (3-13) have uniformly spaced nodes 
and are plotted with a slightly distorted aspect ratio to more readily 
visualize pertinent trends. As previously mentioned these plots represent 
a region of approximately two lamina thicknesses from the free edge. 

Moisture distribution isoclines are numbered on each plot. Iso- 
clines labeled 6 represent fully saturated conditions while those 
labeled 1 represents moisture concentrations of l/6 th the saturated 
values. 


5. RESULTS AND CONCLUSIONS 

The effect of stacking sequence is most noti cable after a short 
Interval, on the order of 2-3 minutes, in a region near the free edge, 
as shown In Figures (3, 5, 8). This phenomena was typical of all 
laminates studied but will be shown using selected laminates. Transverse 
diffusion (z direction) reduces er eliminates the effect of stacking 
sequence for diffusion times greater than five minutes. This was 
indicative of all studies performed, especially in the outer region (see 
Figures (4, 7, 9, 11, 13) of the lamiante, but to a lesser extent in 
the inner region. 

Large transverse gradients in moisture concentrations occurred 
within the outer layers with significant ‘'corner" and "longitudinal" 
gradients in the adjacent and center most layer in all laminates studied, 
(Figure 3-13). Moisture diffusion was more uniform, i.e., smaller 
irregularities, with laminates having lower diffusion rates in the outer 
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most layer. This Is readily seen by comparing the time history contours 
of the [0,±45,90] s and [90,+45,0] s laminates as shown In Figures (10- 
13). 

In conclusion, the stacking sequence was significant for very small 
diffusion times and In a region near the free edge. However, large 
diffusion times tended to eliminate the significance of stacking se- 
quence due to the high transverse diffusion rates and aspect ratio of 
the composite laminate. Moisture concentration gradients within and 
between layers were significant for all diffusion times. From these 
results it is expected that large moisture gradients could produce large 
interlaminar stresses and greatly effect free edge delamination. 
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Figure (2) Longitudinal and Transverse Diffusion 
Coefficients vs. Temperature (°k) 
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Figure (3) Moisture Concentration Contours in a 
[902 >02^5 Laminate at t = 150 sec. 
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Figure (4) Moisture Concentration Contours in a 
[9O 2 0 2 ] $ Laminate at t = 300 sec. 
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Figure (5) Moisture Concentration Contours in a 
[90,0 2 ,90 ] s Laminate at t = 150 sec. 
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Figure (6) Moisture Concentration Contours In a 
[90,02 *90] $ Laminate at t = 200 sec. 
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Figure (7) Moisture Concentration Contours in a 
[90,02 »90] s Laminate at t = 300 sec. 
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Figure (8) Moisture Concentration Contours in a 
[02»902 ] s Laminate at t = 150 sec. 
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Figure (9) Moisture Concentration Contours in a 
[°2»9°2 ^s Laminate at t = 300 sec. 
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Saturated M = 1.4% 
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Figure JO) Moisture Concentration Contours in a 
[0,±45,90] s Laminate at t = 150 sec. 
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Figure (11) Moisture Concentration Contours in a 
[0 t ±45,90] $ laminate at t = 300 sec. 
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Figure (12) Moisture Concentration Contours in a 
[90,±450] $ Laminate at t * 150 sec. 
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Saturated M = 1.4% 
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Figure (13) Moisture Concentration Contours in a 
[90,±45,0] s Laminate at t = 300 sec. 


APPENDIX A 


The finite difference derivatives were developed from a second 
order polynomial with unequal nodal spacing. The coefficients of the 
second order polynomial, (Al), 

2 

y * a 2 x + a^x + a Q (Al) 


can be determined by knowing the magnitude of y at three discrete 
points, x Q , x^, x 2 - Substituting these values of x and y into the 
polynomial, the a^'s can be determined through the matrix equation 
(A2). 



(A2) 


where 


h 5 *1 - x 0 

« = (* 2 - )/(*■) - *„) 


The first derivative at x = x Q in terms of h, a, y Q , y^ , and y 2 is; 


y 0 (° + 2) y^l+a) y 2 
D <V h(l+a) * fia ha(l+a) 


(A3) 


In a similar fashion the first and second central difference deriva- 
tives are; 


D(y 0 ) 



(A4) 
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along z = 0 where 1*1* and 


-"i.ih.i + 2) 


+M i.2 a. , ' “ijd + « i(i ) 


0 + «,,,) 


M 


LI 


= 0 


(A8) 


along y = 0 where j - 1. 

The coefficients used in equations (A6) - (A8) are defined as; 

M. . = moisture concentration at the i,j th node. 

■ *J 

h. . = horizontal step size, (in the y direction) 

* »J 

k. . = vertical step size, (in the z direction) 

* *J 

= diffusion coefficient in the y direction at node i,j 
determined from moisture and temperature data at time t. 


i»j,t 


* D(T,M,Y,Z,e) = D 22 cos 2 0 + sin 2 0 . 


D 11* D 22 5 longitudinal and transverse differsion constants relative 
to the principal material coordinates. 

^ L 

D s diffusion constant in the z direction at the i,j node 

Z i , J ,t 

determined from moisture and temperature data at time t. 
In general D z = D 22 . 

a. . = ratio of distance between i-1, j-i and i, j-i+1, j nodes. 

= ratio of distance between i, y-l-i, j and i, j-i, j+! 
nodes. 

= time increment 


itJ 




APPENDIX B 


The unmodified Gauss Siedel simultaneous equation algorithm 
to solve equation (7) is of the form 


{M k+1 > = (D) - [B]{M k } 


(Bl) 


where 


D i * V c i,i 


B. . = C. ./C- • for i f j 

i,j i,j i,i J 

B. - = 0, with i , j = 1,2,... number of unknowns (NNUM) 

■ t • 


and {M k > is continuously updated from the computed results of {M k+ j}. 

This algorithm, as posed, requires storage on the order of 
(NNUM) and multiplications equal to (NNUM) , per convergence iteration. 
Upon formulating the diffusion problem in conjunction with the Gauss 
Siedel solution scheme there are, at most, four non-zero B- . coef- 

' »J 

ficients per row. The number of non-zero coefficients is a function 
of the order of the difference equation and its order of accuracy. 

By eliminating the large number of zero coefficients in [B] and their 
multiplications the computer resources used to solve this problem 
would be greatly reduced. This was accomplished by developing only 
the non-zero [B] coefficients along with a series of pointer vectors 
relating them to the appropriate guess solution in vector {M^} . 

As coded in HYDIP coefficient storage, including pointer vectors, is 


25 


26 


9*(NNUM) locations while multiplications per convergence iteration 
equals 4*(NNUM). 

Table (Bl) presents a comparison of required computer resources 
to solve the diffusion problem using the modified Gauss Siedel, un- 
modified Gauss Siedel and a non-banded or blocked out of core 
Gaussian elimination algorithms. As seen from the table, the re- 
quired computer resources varies linearly with respect to the number 
of unknowns for the modified Gauss Siedel algorithm while the other 
methods are quadratic or cubic functions of the number of unknowns. 


Table Bl. Storage and Computation Comparison of Three Solution 
Algorithms. 



Number of 
Unknowns 

Number of Mult, 
and Divisions 

Number of Storage 
Locations 

Modified 
Gauss Siedel 


400* 

900 

500 

2000* 

4500 

1000 

4000* 

9000 

Gauss Siedel 

100 

1 x 10 6 * 

10000 

500 

1.25 x 10 8 * 

2.5 x 10 5 

1000 


1.0 x 10 6 

Gaussian 

Elimination 

o 

o 

2 x 10 5 

10000 

500 

1.42 x 10 6 

2.5 x 10 5 

1000 



8.9 x 10 7 

1.0 x 10 6 


* for the Gauss Siedel and the modified Gauss Siedel routines this 
represents the number of multiplications per convergence 
iteration. 






















